;-----------------------------------------------------------------
; 
; bivariate PDF distribution of RN222, SO2, and T  
;
;   - Using a user-defined function for cleaner code
;   - Subselecting an array using its coordinate arrays
;   - Using the "where" function
;   - Using "raster" contour mode for better delination of data
;   - Using "plt_pdfxy"
;-----------------------------------------------------------------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
load "$KAILIB_ROOT/lib_radon.ncl" 


begin 

  limit_XY = False 
  limit_XY = True 

  Lev = "L30"

  load "load_lev.ncl" 

  StringFontHeightF  = 0.030
  cnLineLabelFontHeightF = 0.025
  tmXBLabelFontHeightF = 0.025
  tmYLLabelFontHeightF = 0.025
  tiMainFontHeightF  = 0.030
  tiXAxisFontHeightF = 0.030
  tiYAxisFontHeightF = 0.030
  lbLabelFontHeightF = 0.025

  gsnStringFont = "helvetica-bold"
  gsnStringFont = "helvetica"

  np = 12

  abcd = (/"a","b","c","d","e","f","g","h","i","j","k","l","m","n"/)

  abcd = abcd + ") "

  datapath = "./" 
  plotpath = "./" 
  plotform = "eps"     

  dist = "RN222_SO2" 
  dist = "RN222_T_P" 
  dist = "RN222_T" 

  sitename = "Freiburg" 
  sitename = "Beijing"
  sitename = "Schauinsland" 
  sitename = "Socorro" 
  sitename = "EUR" 
  sitename = "ASIA2" 
  sitename = "USA" 
  sitename = "AUS" 
  sitename = "PANNEL"+"_"+Lev 

  plotname = dist + "_" + sitename 

  flnma = "PDF_CHINA"+Lev+".nc" 
  flnmb = "PDF_EUR"+Lev+".nc"
  flnmc = "PDF_USA"+Lev+".nc" 
  flnmd = "PDF_RUS"+Lev+".nc" 
  flnmo = "PDF_pannel.nc" 


  fla = addfile(datapath+flnma, "r")
  flb = addfile(datapath+flnmb, "r")
  flc = addfile(datapath+flnmc, "r")
  fld = addfile(datapath+flnmd, "r")

  pdfaa = fla->PDFA
  pdfab = fla->PDFB
  pdfac = fla->PDFC

  pdfba = flb->PDFA
  pdfbb = flb->PDFB
  pdfbc = flb->PDFC

  pdfca = flc->PDFA
  pdfcb = flc->PDFB
  pdfcc = flc->PDFC


  pdfda = fld->PDFA
  pdfdb = fld->PDFB
  pdfdc = fld->PDFC


  system("rm "+flnmo) 

  flo = addfile(flnmo,"c") 

  flo->PDFAA=pdfaa 
  flo->PDFAB=pdfab 
  flo->PDFAC=pdfac

  flo->PDFBA=pdfba 
  flo->PDFBB=pdfbb 
  flo->PDFBC=pdfbc

  flo->PDFCA=pdfca 
  flo->PDFCB=pdfcb 
  flo->PDFCC=pdfcc

  flo->PDFDA=pdfda
  flo->PDFDB=pdfdb
  flo->PDFDC=pdfdc

;-----------------------------------------------------------------
; plots
;-----------------------------------------------------------------

  colormap = "BlueDarkRed18" 
  colormap = "BlueDarkOrange18" 
  colormap = "gui_default" 
  colormap = "amwg" 
  colormap = "amwg_blueyellowred" 
  colormap = "WhViBlGrYeOrRe" 
  colormap = "testcmap" 

  plot = new(9,"graphic")

  wks = gsn_open_wks(plotform, plotname)

  gsn_define_colormap(wks,colormap)

  res                      = True
  res@gsnDraw              = False    ; do not draw picture
  res@gsnFrame             = False    ; do not advance frame
  res@cnFillOn             = True     ; turn on color fill
  res@cnFillMode           = "RasterFill"       ; Raster Mode
  res@cnLinesOn            = False    ; no contour lines
  res@cnLineLabelsOn       = False    ; no contour line labels
  res@cnInfoLabelOn        = False
  res@gsnSpreadColors      = True     ; use full colormap

  if(colormap.eq."WhViBlGrYeOrRe") then 
    res@gsnSpreadColorStart  = 15                    ; start at color 10
    res@gsnSpreadColorEnd    = 80                    ; end at color 96
  end if 

  if(colormap.eq."testcmap") then 
    res@gsnSpreadColorStart  = 48                    ; start at color 10
    res@gsnSpreadColorEnd    = 145                    ; end at color 96
    res@gsnSpreadColorEnd    = 175                    ; end at color 96
  end if 

 ;res@cnLabelBarEndStyle   = "ExcludeOuterBoxes"

 ;res@lbLabelBarOn         = False
  res@lbLabelAutoStride    = True
  res@lbLabelAutoStride    = True
  res@lbOrientation        = "vertical" 
  res@lbLabelFontHeightF = lbLabelFontHeightF

  res@tiXAxisString        = "IPRR (cm~S~-3~N~ s~S~-1~N~ ) " 
  res@tiYAxisString        = "Temperature (K) " 
  ;res@tiXAxisOffsetYF      = 0.072 


  if(limit_XY) then 
  res@trXMaxF              = 20.
  res@trYMaxF              = 290.
  res@trYMinF              = 250.
  end if 

;;----------------------------------------------------------------------
;; strings
;;----------------------------------------------------------------------
  res@gsnLeftString    = ""
  res@gsnCenterString  = ""
  res@gsnRightString   = ""

  res@gsnLeftStringFontHeightF   = StringFontHeightF
  res@gsnCenterStringFontHeightF = StringFontHeightF
  res@gsnRightStringFontHeightF  = StringFontHeightF

  res@gsnStringFont = gsnStringFont

;;----------------------------------------------------------------------
;; tickmark
;;----------------------------------------------------------------------
;;res@tmYLMode = "Explicit"
;;res@tmYLValues = (/100.,200.,400.,600.,800.,1000.,2000.,3000.,4000./)
;;res@tmYLLabels = (/"100","200","400","600","800","1000","2000","3000","4000"/)
;;res@tmYLMinorValues  = ispan(0,4000,200)

  res@tmXBLabelFontHeightF = tmXBLabelFontHeightF
  res@tmYLLabelFontHeightF = tmYLLabelFontHeightF

;;----------------------------------------------------------------------
;; axis
;;----------------------------------------------------------------------
  res@tiXAxisFontHeightF = tiXAxisFontHeightF
  res@tiYAxisFontHeightF = tiYAxisFontHeightF


  res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       =  0.2               ; set min contour level
  res@cnMaxLevelValF       =  2.                ; set max contour level
  res@cnLevelSpacingF      =  0.2               ; set contour spacing
  ;res@cnMaxLevelValF       =  1.                ; set max contour level
  ;res@cnLevelSpacingF      =  0.1               ; set contour spacing

  res@gsnLeftString        = abcd(0) + "China, merged"
  res@gsnCenterString      = "" 
  res@gsnRightString       = "DJF" 
  plot(0) = gsn_csm_contour (wks, pdfac, res)
 ;plot(0) = gsn_contour (wks, pdfac, res)

  res@gsnLeftString        = abcd(1) + "USA, merged"
  res@gsnCenterString      = ""
  res@tiYAxisString        = "" 
  plot(1) = gsn_csm_contour (wks, pdfcc, res)
 ;plot(1) = gsn_contour (wks, pdfcc, res)

  res@gsnLeftString        = abcd(2) + "Europe, merged"
  res@gsnCenterString      = ""
  res@tiYAxisString        = "" 
  plot(2) = gsn_csm_contour (wks, pdfbc, res)
 ;plot(2) = gsn_contour (wks, pdfbc, res)

  if(limit_XY) then 
  res@trYMaxF              = 270.
  res@trYMinF              = 220.
  end if 

  res@gsnLeftString        = abcd(3) + "Russia, SW1998 scaled"
  res@gsnCenterString      = ""
  res@tiYAxisString        = "Temperature (K) " 
  plot(3) = gsn_csm_contour (wks, pdfdb, res)
 ;plot(3) = gsn_contour (wks, pdfdb, res)

  res@gsnLeftString        = abcd(4) + "Russia, WCRP1995"
  res@gsnCenterString      = ""
  res@tiYAxisString        = "" 
  plot(4) = gsn_csm_contour (wks, pdfda, res) 
 ;plot(4) = gsn_contour (wks, pdfda, res) 

  res@gsnLeftString        = abcd(5) + "Russia, merged"
  res@gsnCenterString      = ""
  res@tiYAxisString        = "" 
  plot(5) = gsn_csm_contour (wks, pdfdc, res) 
 ;plot(5) = gsn_contour (wks, pdfdc, res) 


  resP                     = True                ; modify the panel plot
 ;resP@gsnMaximize         = True        ; must include w/ Paper Orientation
 ;resP@txString            = title
 ;resP@gsnPanelRowSpec     = True                   ; tell panel what order to plt
 ;resP@gsnPanelLabelBar    = True
 ;resP@lbLabelAutoStride   = True
 ;resP@lbOrientation       = "vertical"          ; vertical label bar
  resP@gsnPanelYWhiteSpacePercent = 5
 ;resp@gsnPanelXWhiteSpacePercent = 5


 ;gsn_panel(wks,plot,(/3,3/),resP)
  gsn_panel(wks,plot,(/3,3/),resP)

;-----------------------------------------------------------------
; Use "plt_pdfxy" to plot
;-----------------------------------------------------------------

;  res@gsnCenterString      = "PDF"
;  plot(0) = plt_pdfxy (wks, pdfa, res)
;
;  resP                     = True 
;  resP@txString            = "" 
;  resP@gsnPanelRowSpec     = True              
;  gsn_panel(wks,plot,(/1,1/),resP)

  ploteps = plotname + ".eps" 
  plotpdf = plotname + ".pdf" 

  system("ps2pdf "+ploteps) 

end                 

